using Plots
using DifferentialEquations
using LaTeXStrings
using Roots
using Graphs
using StatsBase
using GraphPlot
using GraphRecipes
The SIR model is described by the set of equations
$$ \frac{dS(t)}{dt} = - \beta S(t)I(t)$$
$$ \frac{dI}{dt} = \beta S(t)I(t) - rI(t) $$
$$ \frac{dR}{dt} = rI(t) $$
Where
$$ R_0 = \frac{\beta N}{r}$$
function SIR(dudt , u , Param , t)
S,I,R = u
β , r = Param
dudt[1] = - β * S * I
dudt[2] = β * S * I - r * I
dudt[3] = r * I
end
function R0(β,N,r)
return β * N / r
end
function solve_SIR(Param, u0, Tf, nit::Int)
S = zeros(nit)
I = zeros(nit)
R = zeros(nit)
t = range(0, Tf, length=nit)
tspan = (0.0, Tf)
prob = ODEProblem(SIR, u0, tspan, Param)
sol = solve(prob, saveat=t)
for q in 1:length(sol.u)
S[q] = sol.u[q][1]
I[q] = sol.u[q][2]
R[q] = sol.u[q][3]
end
return S, I, R, t
end
solve_SIR (generic function with 1 method)
N = 1000
Tf = 100
Param = [1/N,.1]
tspan = (0.0, Tf)
u0 = [N-1,1,0]
prob = ODEProblem(SIR , u0 , tspan , Param)
sol = solve(prob,tstops=0.0:0.1:Tf)
S = [ q[1] for q in sol.u]
I = [ q[2] for q in sol.u]
R = [ q[3] for q in sol.u]
t = range(0,Tf , length(S))
plot(t , S ./ N , label = L"\frac{S(t)}{N}",size=(1200,800) )
plot!(t , I ./ N, label = L"\frac{I(t)}{N}")
plot!(t , R ./ N, label = L"\frac{R(t)}{N}")
xlabel!(L"t")
title!("β = $(1/N) , r = $(.1) , R0 = $(R0(1/N,N,.1))" , label =L"S(t)")
N = 100
Tf = 1000.0
nit = 1000
tspan = (0.0, Tf)
ϵ = .0001
R_ = [.05,.1,.15,.2]
B = [q/N - ϵ for q in R_]
r0 = [R0(B[q],N,R_[q]) for q in 1:length(R_)]
S = zeros(nit,length(R_))
I = zeros(nit,length(R_))
R = zeros(nit,length(R_))
t = range(0,Tf,nit)
for q in 1:length(R_)
β = B[q]
r = R_[q]
Param = [β,r]
u0 = [N-1,1,0]
prob = ODEProblem(SIR , u0 , tspan , Param)
sol = solve(prob,saveat=range(0,Tf,nit))
for i in 1:length(sol.u)
S[i,q] = sol.u[i][1]
I[i,q] = sol.u[i][2]
R[i,q] = sol.u[i][3]
end
end
gr()
plot1 = plot(t,S[1:end,1] , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])" , label =L"S(t)")
plot!(t,I[1:end,1] , label = L"I(t)")
plot!(t,R[1:end,1] , label = L"R(t)")
plot2 = plot(t,S[1:end,2] , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])", label =L"S(t)")
plot!(t,I[1:end,2], label = L"I(t)")
plot!(t,R[1:end,2], label = L"R(t)")
plot3 = plot(t,S[1:end,3] , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])", label =L"S(t)")
plot!(t,I[1:end,3], label = L"I(t)")
plot!(t,R[1:end,3], label = L"R(t)")
plot4 = plot(t,S[1:end,4] , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])", label =L"S(t)")
plot!(t,I[1:end,4], label = L"I(t)")
plot!(t,R[1:end,4], label = L"R(t)")
plot(plot1,plot2,plot3,plot4 ,size=(1200,800) , xlabel = L"t")
Phase Portrait¶
β = B[1]
r = R_[1]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig1 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[2]
r = R_[2]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig2 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[3]
r = R_[3]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig3 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[4]
r = R_[4]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig4 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)" , label = false , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])")
for j in 10:10:N
plot!(I,S , label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
plot(fig1,fig2,fig3,fig4, size = (1200,1000))
Critical Case $R_0$ = 1¶
N = 100
Tf = 1000.0
nit = 1000
tspan = (0.0, Tf)
ϵ = 0
R_ = [.05,.1,.15,.2]
B = [q/N - ϵ for q in R_]
r0 = [R0(B[q],N,R_[q]) for q in 1:length(R_)]
S = zeros(nit,length(R_))
I = zeros(nit,length(R_))
R = zeros(nit,length(R_))
t = range(0,Tf,nit)
for q in 1:length(R_)
β = B[q]
r = R_[q]
Param = [β,r]
u0 = [N-1,1,0]
prob = ODEProblem(SIR , u0 , tspan , Param)
sol = solve(prob,saveat=range(0,Tf,nit))
for i in 1:length(sol.u)
S[i,q] = sol.u[i][1]
I[i,q] = sol.u[i][2]
R[i,q] = sol.u[i][3]
end
end
plot1 = plot(t,S[1:end,1] , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])" , label =L"S(t)")
plot!(t,I[1:end,1], label = L"I(t)")
plot!(t,R[1:end,1], label = L"R(t)")
plot2 = plot(t,S[1:end,2] , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])", label =L"S(t)")
plot!(t,I[1:end,2], label = L"I(t)")
plot!(t,R[1:end,2], label = L"R(t)")
plot3 = plot(t,S[1:end,3] , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])", label =L"S(t)")
plot!(t,I[1:end,3], label = L"I(t)")
plot!(t,R[1:end,3], label = L"R(t)")
plot4 = plot(t,S[1:end,4] , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])", label =L"S(t)")
plot!(t,I[1:end,4], label = L"I(t)")
plot!(t,R[1:end,4], label = L"R(t)")
plot(plot1,plot2,plot3,plot4 ,size=(1200,800) , xlabel = L"t")
Phase Portrait¶
β = B[1]
r = R_[1]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig1 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[2]
r = R_[2]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig2 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[3]
r = R_[3]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig3 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[4]
r = R_[4]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig4 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)" , label = false , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])")
for j in 10:10:N
plot!(I,S , label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
plot(fig1,fig2,fig3,fig4, size = (1200,1000))
Supercritical case $R_0 > 1$¶
N = 100
Tf = 1000.0
nit = 9000
tspan = (0.0, Tf)
ϵ = 0.0008
R_ = [.05,.1,.15,.2]
B = [q/N + ϵ for q in R_]
r0 = [R0(B[q],N,R_[q]) for q in 1:length(R_)]
S = zeros(nit,length(R_))
I = zeros(nit,length(R_))
R = zeros(nit,length(R_))
t = range(0,Tf,nit)
for q in 1:length(R_)
β = B[q]
r = R_[q]
Param = [β,r]
u0 = [N-1,1,0]
prob = ODEProblem(SIR , u0 , tspan , Param)
sol = solve(prob,saveat=range(0,Tf,nit))
for i in 1:length(sol.u)
S[i,q] = sol.u[i][1]
I[i,q] = sol.u[i][2]
R[i,q] = sol.u[i][3]
end
end
plot1 = plot(t,S[1:end,1] , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])" , label =L"S(t)")
plot!(t,I[1:end,1], label = L"I(t)")
plot!(t,R[1:end,1], label = L"R(t)")
plot2 = plot(t,S[1:end,2] , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])", label =L"S(t)")
plot!(t,I[1:end,2], label = L"I(t)")
plot!(t,R[1:end,2], label = L"R(t)")
plot3 = plot(t,S[1:end,3] , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])", label =L"S(t)")
plot!(t,I[1:end,3], label = L"I(t)")
plot!(t,R[1:end,3], label = L"R(t)")
plot4 = plot(t,S[1:end,4] , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])", label =L"S(t)")
plot!(t,I[1:end,4], label = L"I(t)")
plot!(t,R[1:end,4], label = L"R(t)")
plot(plot1,plot2,plot3,plot4 ,size=(1200,800) , xlabel = L"t")
Phase Portrait¶
β = B[1]
r = R_[1]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig1 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[1]) , r = $(R_[1]) , R0 = $(r0[1])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[2]
r = R_[2]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig2 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[2]) , r = $(R_[2]) , R0 = $(r0[2])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[3]
r = R_[3]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig3 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)", label = false , title = "β = $(B[3]) , r = $(R_[3]) , R0 = $(r0[3])")
for j in 10:10:N
plot!(I,S, label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
β = B[4]
r = R_[4]
Param = [β,r]
u0 = [N-1,1,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
fig4 = plot(I,S , ylabel = L"S(t)" , xlabel = L"I(t)" , label = false , title = "β = $(B[4]) , r = $(R_[4]) , R0 = $(r0[4])")
for j in 10:10:N
plot!(I,S , label = false)
u0 = [N-j,j,0]
S , I , R = solve_SIR(Param, u0, Tf, nit::Int)
end
plot!( range(0,N) , N .- range(0,N) , label = L"S(0) = N - I(0)")
plot(fig1,fig2,fig3,fig4, size = (1200,1000))
We can write solution in an implicit form for $S(t) , I(t) , R(t)$
$$ S(t) = S(0) exp \left( -\int_0^t \beta I(s) ds \right)$$
$$ I(t) = I(0) exp \left( \int_0^t \beta (I(s) - r) ds\right) $$
$$ R(t) = R(0) + \int_0^t r I(s) ds $$
$$ S(t) + I(t) + R(t) = N$$
We notice that
$$ S(t) = S(0) exp\left( \frac{-\beta}{r} (R(t) - R(0))\right)$$
We can write
$$ \frac{dR}{dt} = r (N - S(t) - R(t)) $$
Then
$$ \frac{dR(t)}{dt} = r (N - R(t) - S(0) exp\left( \frac{-\beta}{r} (R(t) - R(0))\right) ) $$
We have $ \frac{dR(t)}{dt} \rightarrow 0 $ as $t \rightarrow \infty$. Then we have the following fix-point equation
$$ R(t) = N - S(0) exp\left( \frac{-\beta}{r} (R(t) - R(0))\right) $$
The initial conditions are
$$ \frac{S(0)}{N} = 1 - \epsilon $$
$$ \frac{I(0)}{N} = \epsilon $$
$$ R(0) = 0$$
Then we have
$$ r(t) = 1 - (1-\epsilon) e^{-R_0 r} $$
$$ \frac{dr}{dR_0} = \frac{(1-\epsilon)re^{-R_0 r}}{1 - (1-\epsilon)R_0e^{-R_0r}}$$
r = range(0,1,100)
R0_val = range(0.5 , 2.0 , 4)
R_0 = range(.0,10,1000)
Ε = [.001,.005,.010,.015]
function rt(r,ϵ,R0)
return 1 - (1-ϵ) * exp( - R0 * r)
end
function r_fixed_point(ϵ,r_0)
return find_zero(x -> rt(x,ϵ , r_0) - x, .5)
end
r1 = r_fixed_point( Ε[1] , R0_val[1] )
r2 = r_fixed_point( Ε[2] , R0_val[2] )
r3 = r_fixed_point( Ε[3] , R0_val[3] )
r4 = r_fixed_point( Ε[4] , R0_val[4] )
r_fixed = r_fixed_point.(.001 , R_0)
plot1 = plot(r,r , label = L"y=r(t)" , xlabel = L"r(t)" , size = (1200,800))
plot!(r , rt.(r, Ε[1] ,R0_val[1]) , label = "ϵ = $(Ε[1]) , R0 = $(R0_val[1]) , r* = $(r1)")
plot!(r , rt.(r, Ε[2] ,R0_val[2]) , label = "ϵ = $(Ε[2]) , R0 = $(R0_val[2]) , r* = $(r2)")
plot!(r , rt.(r, Ε[3] ,R0_val[3]) , label = "ϵ = $(Ε[3]) , R0 = $(R0_val[3]) , r* = $(r3)")
plot!(r , rt.(r, Ε[4] ,R0_val[4]) , label = "ϵ = $(Ε[4]) , R0 = $(R0_val[4]) , r* = $(r4)")
scatter!([r1],[r1] , label = false)
scatter!([r2],[r2] , label = false)
scatter!([r3],[r3] , label = false)
scatter!([r4],[r4] , label = false)
plot2 = plot(R_0 , r_fixed , xlabel = L"R_0" , title = "ϵ = $(0.001)" , label = L"r^{*}")
scatter!([1] , [r_fixed_point( .001 , 1.0 )] , label = "R0 = 1 , r* = $(r2)")
plot(plot1,plot2)
drdR0 = (r_fixed[2:end] - r_fixed[1:end-1]) ./ (R_0[2:end] - R_0[1:end-1])
plot(R_0[1:end-1] , drdR0 , size = (1200,800) , title = L"\frac{dr^{*}}{dR_0}" , xlabel = L"R_0")
Exercice II¶
function Infection( G::Graph , p , nit , I0=1)
### Infection with Recovery
nodes = nv(G)
S_Matrix = ones(Int , nit , nv(G) )
I_Matrix = zeros(Int , nit , nv(G) )
R_Matrix = zeros(Int , nit , nv(G) )
S_Matrix[1:end,I0] .= 0
I_Matrix[1,I0] = 1
R_Matrix[2:end,I0] .= 1
for q in 1:nit-1
for j in 1:nodes
if (I_Matrix[q,j] == 1) && (R_Matrix[q,j] == 0)
for k in neighbors(G , j)
if( rand() < p ) && (R_Matrix[q,k] == 0) && (I_Matrix[q,k] == 0)
S_Matrix[q+1:end,k] .= 0
I_Matrix[q+1,k] = 1
R_Matrix[q+2:end,k] .= 1
end
end
end
end
end
return S_Matrix , I_Matrix , R_Matrix
end
function Infection_sim(G::Graph , p , nit , nsim,I0=1)
S_av = zeros(nit)
I_av = zeros(nit)
R_av = zeros(nit)
for q in 1:nsim
S , I , R = Infection(G,p,nit,I0)
S_av[1:end] .= S_av[1:end] .+ sum(S,dims=2) / nsim
I_av[1:end] .= I_av[1:end] .+ sum(I,dims=2) / nsim
R_av[1:end] .= R_av[1:end] .+ sum(R,dims=2) / nsim
end
return S_av , I_av , R_av
end
function statistics( G::Graph , p , nit , nsim)
N = nv(G)
I_max = 0
prop_infect = 0
clear_I = 0
nodes = [q for q in 1:N]
for q in 1:nsim
I0 = rand(nodes)
S_matrix , I_matrix , R_matrix = Infection(G,p,nit,I0)
S = vec( sum(S_matrix , dims = 2) )
I = vec( sum(I_matrix , dims = 2) )
R = vec( sum(R_matrix , dims = 2) )
clear_I += findfirst(x -> x ==0 , I) / nsim
I_max += argmax(I) / nsim
prop_infect += (R[end] / N) / nsim
end
return [ prop_infect , I_max , clear_I]
end
statistics (generic function with 1 method)
N = 100 ### Number of nodes
nit = 200 ### Number of Iterations
nsim = 500 ### Number of Simulations
p1 = .4 ### Probability of Infection
p2 = .45 ### Probability of Infection
p3 = .5 ### Probability of Infection
p4 = .55 ### Probability of Infection
0.55
Square Graph¶
G1 = Graphs.grid([Int(sqrt(N)) , Int(sqrt(N))])
S1 , I1 , R1 = Infection_sim(G1 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim(G1 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim(G1 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim(G1 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) , xlabel = "Iteration")
Random Graph¶
G2 = watts_strogatz(N , 4 , 1)
S1 , I1 , R1 = Infection_sim(G2 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim(G2 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim(G2 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim(G2 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) , xlabel = "Iteration")
Watts Strogaz Network¶
G3 = watts_strogatz(N , 4 , .1)
S1 , I1 , R1 = Infection_sim(G3 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim(G3 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim(G3 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim(G3 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) , xlabel = "Iteration")
Barabasi Albert model¶
G4 = barabasi_albert(N , 2)
S1 , I1 , R1 = Infection_sim(G4 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim(G4 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim(G4 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim(G4 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) , xlabel = "Iteration" )
c) We see that the discrete model reproduces welll the continuous curves. However, the model becomes strongly dependent on the network topological parameters and the rate parameters $\beta,\gamma,R_0$ lose their meaning
d)
Square Graph¶
N = 100 ### Number of nodes
p = range(0,1,100) ### Proabilities Range
nit = 200 ### Number of Iterations
nsim = 5000 ### Number of SImulations
G1 = Graphs.grid([Int(sqrt(N)) , Int(sqrt(N))])
Stats = map( x -> statistics(G1 , x , nit , nsim) , p )
prop_infect = [q[1] for q in Stats]
I_max = [q[2] for q in Stats]
clear_I = [q[3] for q in Stats]
plot1 = plot(p,prop_infect , title = "Average Total Infected Proportion" , label =false)
plot2 = plot(p,I_max , title = "Average Maximum Infected Iteration", label =false)
plot3 = plot(p,clear_I , title = "Average Clear Infection Iteration", label =false)
plot(plot1 , plot2 , plot3, xlabel ="Probability of Infection" , size = (1200,600) , layout=(1,3))
Random Graph¶
N = 200 ### Number of nodes
p = range(0,1,100) ### Proabilities Range
nit = 200 ### Number of Iterations
nsim = 5000 ### Number of SImulations
G2 = watts_strogatz(N , 5 , 1)
Stats = map( x -> statistics(G2 , x , nit , nsim) , p )
prop_infect = [q[1] for q in Stats]
I_max = [q[2] for q in Stats]
clear_I = [q[3] for q in Stats]
plot1 = plot(p,prop_infect , title = "Average Total Infected Proportion" , label =false)
plot2 = plot(p,I_max , title = "Average Maximum Infected Iteration", label =false)
plot3 = plot(p,clear_I , title = "Average Clear Infection Iteration", label =false)
plot(plot1 , plot2 , plot3, xlabel ="Probability of Infection" , size = (1200,600) , layout=(1,3))
Watts-Strogatz Graph¶
N = 200 ### Number of nodes
p = range(0,1,100) ### Proabilities Range
nit = 200 ### Number of Iterations
nsim = 5000 ### Number of SImulations
G3 = watts_strogatz(N , 5 , .1)
Stats = map( x -> statistics(G3 , x , nit , nsim) , p )
prop_infect = [q[1] for q in Stats]
I_max = [q[2] for q in Stats]
clear_I = [q[3] for q in Stats]
plot1 = plot(p,prop_infect , title = "Average Total Infected Proportion" , label =false)
plot2 = plot(p,I_max , title = "Average Maximum Infected Iteration", label =false)
plot3 = plot(p,clear_I , title = "Average Clear Infection Iteration", label =false)
plot(plot1 , plot2 , plot3, xlabel ="Probability of Infection" , size = (1200,600) , layout=(1,3))
Albert-Barabasi Graph¶
N = 200 ### Number of nodes
p = range(0,1,100) ### Proabilities Range
nit = 200 ### Number of Iterations
nsim = 5000 ### Number of SImulations
G4 = barabasi_albert(N , 2)
Stats = map( x -> statistics(G4 , x , nit , nsim) , p )
prop_infect = [q[1] for q in Stats]
I_max = [q[2] for q in Stats]
clear_I = [q[3] for q in Stats]
plot1 = plot(p,prop_infect , title = "Average Total Infected Proportion" , label =false)
plot2 = plot(p,I_max , title = "Average Maximum Infected Iteration", label =false)
plot3 = plot(p,clear_I , title = "Average Clear Infection Iteration", label =false)
plot(plot1 , plot2 , plot3, xlabel ="Probability of Infection" , size = (1200,600) , layout=(1,3))
e)
f)
N = 25 ### Number of nodes
nit = 25 ### Number of Iterations
p = .7 ### Infection Probability
0.7
Square Graph¶
G1 = Graphs.grid([Int(sqrt(N)) , Int(sqrt(N))])
S , I , R = Infection(G1 , p , nit)
Colors = [ R[q,j] == 1 ? :red : :blue for q in 1:nit , j in 1:nv(G1)]
plot(G1, nodecolor=Colors[1,1:end])
anim = @animate for q in 2:nit
plot(G1, nodecolor=nodecolor=Colors[q,1:end],title="Itreation = $(q)")
end
gif(anim, fps=1)
┌ Info: Saved animation to /tmp/jl_9wvYbmlGty.gif └ @ Plots /home/lucas/.julia/packages/Plots/a3u1v/src/animation.jl:156
Random Graph¶
G2 = watts_strogatz(N , 4 , 1)
S , I , R = Infection(G2 , p , nit)
Colors = [ R[q,j] == 1 ? :red : :blue for q in 1:nit , j in 1:nv(G1)]
plot(G2, nodecolor=Colors[1,1:end])
anim = @animate for q in 2:nit
plot(G2, nodecolor=nodecolor=Colors[q,1:end],title="Itreation = $(q)")
end
gif(anim, fps=1)
┌ Info: Saved animation to /tmp/jl_nArM1i3k7O.gif └ @ Plots /home/lucas/.julia/packages/Plots/a3u1v/src/animation.jl:156
Watts-Strogatz¶
G3 = watts_strogatz(N , 4 , .1)
S , I , R = Infection(G1 , p , nit)
Colors = [ R[q,j] == 1 ? :red : :lightblue for q in 1:nit , j in 1:nv(G1)]
plot(G3, nodecolor=Colors[1,1:end])
anim = @animate for q in 2:nit
plot(G3, nodecolor=nodecolor=Colors[q,1:end],title="Itreation = $(q)")
end
gif(anim, fps=1)
┌ Info: Saved animation to /tmp/jl_1aNy969YZU.gif └ @ Plots /home/lucas/.julia/packages/Plots/a3u1v/src/animation.jl:156
Albert-Barabasi¶
G4 = barabasi_albert(N , 2)
S , I , R = Infection(G1 , p , nit)
Colors = [ R[q,j] == 1 ? :red : :lightblue for q in 1:nit , j in 1:nv(G1)]
plot(G4, nodecolor=Colors[1,1:end])
anim = @animate for q in 2:nit
plot(G4, nodecolor=nodecolor=Colors[q,1:end],title="Itreation = $(q)")
end
gif(anim, fps=1)
┌ Info: Saved animation to /tmp/jl_7w36tarkN8.gif └ @ Plots /home/lucas/.julia/packages/Plots/a3u1v/src/animation.jl:156
Exercice III¶
function Infection_asynchronous( G::Graph , p , nit , I0=1)
### Infection with Recovery
nodes = nv(G)
S_Matrix = ones(Int , nit , nv(G) )
I_Matrix = zeros(Int , nit , nv(G) )
R_Matrix = zeros(Int , nit , nv(G) )
S_Matrix[1:end,I0] .= 0
I_Matrix[1,I0] = 1
R_Matrix[2:end,I0] .= 1
for q in 1:nit-1
I = rand(1:nodes)
if( I_Matrix[q,I] == 1)
for k in neighbors(G,I)
if( rand() < p ) && (R_Matrix[q,k] == 0) && (I_Matrix[q,k] == 0)
S_Matrix[q+1:end,k] .= 0
I_Matrix[q+1,k] = 1
R_Matrix[q+2:end,k] .= 1
end
end
else
S_Matrix[q+1,1:end] .= S_Matrix[q,1:end]
I_Matrix[q+1,1:end] .= I_Matrix[q,1:end]
R_Matrix[q+1,1:end] .= R_Matrix[q,1:end]
end
end
return S_Matrix , I_Matrix , R_Matrix
end
function Infection_sim_asynchronous(G::Graph , p , nit , nsim,I0=1)
S_av = zeros(nit)
I_av = zeros(nit)
R_av = zeros(nit)
for q in 1:nsim
S , I , R = Infection_asynchronous(G,p,nit,I0)
S_av[1:end] .= S_av[1:end] .+ sum(S,dims=2) / nsim
I_av[1:end] .= I_av[1:end] .+ sum(I,dims=2) / nsim
R_av[1:end] .= R_av[1:end] .+ sum(R,dims=2) / nsim
end
return S_av , I_av , R_av
end
Infection_sim_asynchronous (generic function with 2 methods)
N = 100 ### Number of nodes
nit = 2000 ### Number of Iterations
nsim = 500 ### Number of Simulations
p1 = .4 ### Probability of Infection
p2 = .45 ### Probability of Infection
p3 = .5 ### Probability of Infection
p4 = .55 ### Probability of Infection
0.55
Square Graph¶
G1 = Graphs.grid([10,10])
S1 , I1 , R1 = Infection_sim_asynchronous(G1 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim_asynchronous(G1 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim_asynchronous(G1 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim_asynchronous(G1 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800))
Random Graph¶
G2 = watts_strogatz(N , 4 , 1)
S1 , I1 , R1 = Infection_sim_asynchronous(G2 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim_asynchronous(G2 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim_asynchronous(G2 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim_asynchronous(G2 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) )
Watts-Strogaz Graph¶
G3 = watts_strogatz(N , 4 , .1)
S1 , I1 , R1 = Infection_sim_asynchronous(G3 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim_asynchronous(G3 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim_asynchronous(G3 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim_asynchronous(G3 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) )
Albert-Barabasi Graph¶
G4 = barabasi_albert(N , 2)
S1 , I1 , R1 = Infection_sim_asynchronous(G4 , p1 , nit , nsim)
S2 , I2 , R2 = Infection_sim_asynchronous(G4 , p2 , nit , nsim)
S3 , I3 , R3 = Infection_sim_asynchronous(G4 , p3 , nit , nsim)
S4 , I4 , R4 = Infection_sim_asynchronous(G4 , p4 , nit , nsim)
plot1 = plot(S1 / N , label = L"S(k)" , title = "p=$(p1)")
plot!(I1 / N , label = L"I(k)")
plot!(R1 / N , label = L"R(k)")
plot2 = plot(S2 / N , label = L"S(k)" , title = "p=$(p2)")
plot!(I2 / N , label = L"I(k)")
plot!(R2 / N , label = L"R(k)")
plot3 = plot(S3 / N , label = L"S(k)" , title = "p=$(p3)")
plot!(I3 / N , label = L"I(k)")
plot!(R3 / N , label = L"R(k)")
plot4 = plot(S4 / N , label = L"S(k)" , title = "p=$(p4)")
plot!(I4 / N , label = L"I(k)")
plot!(R4 / N , label = L"R(k)")
plot(plot1 , plot2 , plot3 , plot4 , size = (1200 , 800) )